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Abstract: An efficient spectral element (SE) with electric potential degrees of freedom 
(DOF) is proposed to investigate the static electromechanical responses of a piezoelectric 
bimorph for its actuator and sensor functions. A sublayer model based on the piece wise 
linear approximation for the electric potential is used to describe the nonlinear distribution 
of electric potential through the thickness of the piezoelectric layers. An equivalent single 
layer (ESL) model based on first-order shear deformation theory (FSDT) is used to 
describe the displacement field. The Legendre orthogonal polynomials of order 5 are used 
in the element interpolation functions. The validity and the capability of the present SE 
model for investigation of global and local responses of the piezoelectric bimorph are 
confirmed by comparing the present solutions with those obtained from coupled 3-D finite 
element (FE) analysis. It is shown that, without introducing any higher-order electric 
potential assumptions, the current method can accurately describe the distribution of the 
electric potential across the thickness even for a rather thick bimorph. It is revealed that the 
effect of electric potential is significant when the bimorph is used as sensor while the effect 
is insignificant when the bimorph is used as actuator, and therefore, the present study may 
provide a better understanding of the nonlinear induced electric potential for bimorph 
sensor and actuator. 

Keywords: spectral element method; piezoelectric bimorph; electric potential; sublayer; 
piecewise linear 
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1. Introduction 

Piezoelectric materials generate electric potentials in response to mechanical stresses, and 
conversely, produce mechanical movements in response to electric potentials. Therefore, piezoelectric 
materials can be used both as actuators and sensors, they transform electrical energy into mechanical 
energy, and vice versa. To achieve practically meaningful actuation and sensing capabilities, a 
piezoelectric bimorph consisting of two piezoelectric layers is widely used [1,2]. A broad range of 
electromechanical applications have been reported, such as electroacoustic transducers [3,4], medical 
devices [5], microcantilever biosensors [6], and atomic force microscope (AFM) cantilevers [7]. However, 
before piezoelectric bimorphs can be utilized in all these applications, it is first necessary to investigate 
both the global responses and the local responses, e.g., the deflection and the distribution of the electric 
potential across the thickness. 

There have been many theories and models developed for analyzing piezoelectric bimorph 
structures with emphasis on approximating the mechanical displacement and electric potential. By 
carrying out exact 3-D analytical solutions for the simply supported piezoelectric plate [8,9], it is 
shown that the distribution of the electric potential across the thickness is nearly quadratic. This 
implies that the assumption of linear distribution of the electric potential across the thickness adopted 
by many numerical models [10,11] cannot address this nonlinear electric potential. Since exact 3-D 
analytical solutions are not available for more general cases of loading and boundary conditions, the 
introduction of the finite element (FE) method is desirable. A considerable amount of literature has 
been published on the FE analysis of piezoelectric smart structures [12-14]. Among these works, the 
simplest and often used model is the equivalent single layer (ESL) model in which the displacement 
and strain functions are assumed to be continuous through the thickness. There are two main kinds of 
theories used for ESL models. One is the classical laminated plate theory (CLPT) [15,16], and the 
other one is the shear deformation theory, which branches out into first-order shear deformation theory 
(FSDT) [17,18] and higher order shear deformation theory (HSDT) [19,20]. The ESL model is simple 
and capable of predicting the global responses of the bimorph, but it does not account for the nonlinear 
distribution of the electric potential across the thickness. To overcome this shortcoming, the FE model 
using the layer- wise theory [21-24] or the sublayer theory [2,25-28] has been recommended. In the 
latter case, the piezoelectric layer is divided into appropriate number of thin sublayers. For each of 
these sublayers, a linear electric potential distribution across the plate thickness is assumed. It is further 
expected that the quadratic distribution of the electric potential across the plate thickness can be 
accurately approached with more sublayers adopted. 

Generally, accurately simulation of the local responses of the piezoelectric bimorph structures 
would inevitably lead to a very dense FE mesh when using the FE method. Hence, conventional FE 
simulation becomes computationally very inefficient. A more efficient method is the spectral element 
(SE) method which combines the geometric flexibility of FE method with the high accuracy of the 
pseudo spectral method. This method was first presented by Patera in the mid 1980s [29]. In fact, the 
SE method and FE method are closely related and built on the same ideas. The main difference 
between them is that SE method uses orthogonal polynomials, such as Legendre and Cheybysev 
polynomials, in the shape functions. The SE method results naturally in diagonal mass matrices which 
is a distinct advantage over traditional FE method especially for transient analysis. Moreover, to have 
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an accurate simulation with the conventional FE method, a mesh with a large number of elements and 
degrees of freedom (DOFs) is inevitably needed. The SE method, in which the polynomial order is 
increased and the mesh size is decreased, can be used to overcome this problem. The SE method has 
been widely applied to many engineering problems related to acoustics, fluid dynamics and 
seismology [30-35]. Recently, the SE method has been extensively used to investigate the wave 
propagation problems for the purpose of damage detection in structures [36,37]. However, according 
to the authors' best knowledge, the SE method has not been previously used for accurately modeling of 
the through-the-thickness electric potentials for piezoelectric bimorphs. 

For the purpose of accurately representing the mechanical displacement and the electric potential, a 
reasonable choice is to use the ESL model for the mechanical variables and the layer-wise theory or 
the sublayer theory for the electric variables. In the present work, we attempt to combine the merits of 
the SE method and the sublayer model. More specifically, the mechanical variables, i.e., the 
displacements, are described based on FSDT. The electrical variables, i.e., the potentials, are described 
using the sublayer model. SE method is then utilized to deduce the governing equations. Legendre 
orthogonal polynomials are adopted in the interpolation function to improve the accuracy. To validate 
the effectiveness and the capability of the present model, numerical simulations for a simply supported 
piezoelectric bimorph with two different load cases, i.e., a uniform pressure load applied to the top 
surface and a uniform potential applied to the top and bottom surfaces, are carried out. The results 
obtained by the present approach are then compared to those coming from the coupled 3-D FE 
simulations using ABAQUS. The comparisons show the good accuracy and efficiency of SE method 
for modeling of the through-the-thickness electric potentials of the piezoelectric bimorph. 

2. Mathematical Formulation 

2.7. Constitutive Relationships, Displacement and Strain 

A piezoelectric bimorph made of two identical PZT-4 piezoelectric layers, which has been 
investigated by Fernandes [1], is considered here. The PZT-4 layer is assumed to behave in a linear 
orthotropic manner with small displacements and strains. As depicted in Figure 1, both piezoelectric 
layers have the same thickness 0.5 h and are poled in the same direction. The x-y plane of the 
coordinate system x-y-z coincides with the middle plane of the bimorph, and the z axis is defined 
normal to the middle plane following the right-hand rule. This work aims to investigate the problem of 
a simply supported piezoelectric bimorph under a uniform pressure load or an applied electric potential 
in the framework of linear theory of piezoelectricity. Assuming the PZT-4 layers work under 
isothermal conditions, the pyroelectric effects and thermomechanical couplings are not taken into 
account. Consequently, a linear constitutive relationship addressing both the direct and converse 
piezoelectric effects is utilized for the analysis of the piezoelectric bimorph, which can be written as: 

o = C£ - e T E 

D = ez + gE (1) 

where a = [a x a y a z r yz r zx r^f and c = £y e z lyz lzx lxy f represent stress vector 
and strain vector, respectively, e = f E E E f , the electric field vector, D = f D D D f , the 

I x y z \ [ x y z J 
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electric displacement vector, c , the elastic coefficient matrix, g , the dielectric coefficient matrix, and 
e , the piezoelectric stress coefficient matrix. 

Figure 1. Geometry of a piezoelectric bimorph. 




An ESL model adopting the FSDT is adopted to describe the mechanical displacement. The 
displacement field of a piezoelectric bimorph based on FSDT takes on the form [17,18]: 



u x,y,z,t — u(x, y, t) + za{x, y, t) 
v x,y,z,t =v 0, y, t) + zfi(x, y, t) 
w x,y,z,t w(x,y,t) 



(2) 



where u ,v ,w denote the displacements of an arbitrary point on the mid plane z = 0 , a and — j3 
denote the rotations of a transverse normal about the y and x axes, respectively. In the FSDT, the 
transverse shear strains are assumed to be constant with respect to the thickness coordinate. The 
constant state of transverse shear strains across the thickness is a gross approximation of the true strain 
field, which is at least quadratic through the thickness. 
We define: 



U 



U V w 



U 



u v w a (3 



(3) 
(4) 



where U is the displacement vector, and U is a generalized displacement vector. Then Equation (2) 
can be written in matrix form as: 



where: 



u = zu 

1 0 0 z 0 

0 1 0 0 z 

0 0 10 0 



(5) 



(6) 



The infinitesimal strain components associated with the displacements are given by: 
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£ = LU 



(7) 



where L is the derivation operator defined as: 
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(8) 



2.2. Approximations for Displacements 



The Legendre polynomials based SE method can be described as follows: the bimorph is firstly 
discretized using a set of non-overlapping rectangular elements, as in the traditional FE method. Each 
rectangular element, denoted by fi e , is then mapped to a reference element, denoted by 
ft mi : £ e [— x 7] e [— l,l] , using an invertible local mapping. The discretization procedure is 

illustrated in Figure 2. 

Subsequently, a set of nodes, denoted by 77 ■ , are defined in the local coordinate system £ — 77 
of the reference element fi ref as roots of the following polynomial expression: 

a-a^(e) = o 



(9) 



where P N is the yV-th order Legendre polynomial. In fact, the nodes are the 2-D Gauss-Lobatto- 

Legendre (GLL) points. In contrast to the classical FE method, the distribution of nodes is irregular, as 
shown in Figure 2. In the current formulation, the 5-th order Legendre polynomial is chosen, hence 
36 nodes can be specified in the reference element f2 ref , as depicted in Figure 2. 

Figure 2. Discretization of a plate and an example of spectral element. 



O potential DOFs 
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mechanical DOFs 
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The 1-D shape functions at the 1-D GLL points £. are defined as [36]: 

-1 a-aw 



MO 



for i = iV + 1 



(10) 
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An important property of these interpolation functions is the discrete orthogonality expressed as: 



(11) 



where <5 denotes the Kronecker delta. The 2-D shape functions are constructed as a tensor product of 



the 1-D ones: 



*iM,v) = MOhM for ij = 1,...,7V + 1 



(12) 



Figure 3 shows two examples of the 2-D shape functions which indicate that each shape function 
has the value 1 at one node and vanish at all other nodes. 

Figure 3. Selected shape functions for a 36-node spectral element, (a) ^ 32 (^r]) ; 
(b) ^(M. 





Coordinates x and y within each fi e may be uniquely related to £ and rj upon the invertible mapping: 

6 6 

(z(&v),y(&v)) = (13) 

i=l j=l 



where x { . and y { . denote the coordinates of x and y , respectively, of the element nodes ^ . The 
generalized displacements u 9 v 9 w 9 a and /? over an reference element fi ref are discretized by the 
2-D shape functions as: 



6 6 



(14) 



where u. 9 v.., w - 9 ex.. and are the nodal values of the generalized displacements. The discrete 

IJ 5 ZJ ' ZJ ' ZJ ^ZJ fe ^ 

element nodal displacement vector is expressed as: 



T T T 
All Al2 A66 



(15) 



where q^. is the displacement vector of the node 77 : 
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% v v % % A;f (16) 
Substituting Equations (14) into Equation (15) yields: 

U = N M q e (17) 
where N is the displacement shape function matrix which can be expressed as: 



N N ••• N 



(18) 



with: 

= ^SxS (19) 
where I 5x5 is a 5 x 5 identity matrix. Substituting Equation (17) into Equation (7) yields: 

8 = B q e ( 2 °) 
where is strain-displacement matrix which can be written as: 

B „ = LN U (21) 

2.3. Approximations for Electric Potential 

For the purpose of accurately modeling the distribution of the electric potential across thickness, 
each layer of the piezoelectric bimorph is subdivided mathematically into n thinner sublayers. As 
shown in Figure 4, the sublayers are numbered in top-to-bottom order. The z coordinates of the top 
and bottom surfaces of the i-th sublayer are denoted by z i and z { _ x , respectively. In each sublayer, the 
distribution of the electric potential (j> % {z) is assumed to be linear across the thickness such that: 

<j>\z) = N;O l (22) 

where is the interpolation function and 4>* is a column matrix composed of the electric potentials 
at the top and the bottom surfaces of the i — th sublayer, which can be expressed as: 

1 



0 h 



z- — z z — Z- , 

i i—L 



h = Z ; — Z- , , Z. ■ i < Z < Z- 



v z t _ x ^z^z, ( 23) 



=[4-i M (24) 

In this way, the assumption of linear distribution of electric potential across the thickness is used 
not in the whole piezoelectric layer, but in each sublayer instead. As a result, the electric potential is 
approximated as piecewise linear across the thickness and it is expected that the quadratic distribution 
of the electric potential across the bimorph thickness can be approached with more sublayers adopted. 
As mentioned before, the mechanical displacement field is approximated using ESL model based on 
FSDT and the piezoelectric bimorph is discretized using 2-D mesh. To keep the compatibility, 
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each sublayer of the piezoelectric layer is also discretized using the same mesh. Consequently, an 
element potential vector <D e is then introduced in the spectral plate finite element, which is defined as: 

* e =K A - Kf (25) 

The surface potential of the sublayer, 0. , is assumed to be constant over the element and 
<f> Q , <f) v • • * (f) 2n are treated as elemental DOFs, as illustrated in Figure 2. Furthermore, the top and bottom 
surfaces of the piezoelectric layers are always coated with metallic coatings of zero thickness and the 
potentials on the electrodes should be taken as independent of x, y . Thus the present method combines 
an ESL theory for the displacement and a piecewise linear approximation for the electric potential. 
Under the quasi-electrostatic approximation, the electric field and the electric potential in each 
sublayer have the following relationship: 

W(z) = -B;0* (26) 
where E l (z) is the electric field of the i-th sublayer, is the electric field-potential matrix, given by: 

Bl = VN; (27) 



Figure 4. A sublayer model for a piezoelectric bimorph. 
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2.4. Governing Equations 

By applying Hamilton's principle, the elementary dynamic equations for the piezoelectric bimorph 
plate can be obtained: 



M^q 6 + K e „„q e + K^O e = F, 
Klq e + K^O e = F 



(28) 



where M e denotes the element mass matrix; K e , mechanical stiffness matrix; Kf , and K e , the 

uu ' uu ' ' u<p (pu 



piezoelectric coupling matrices; the dielectric permittivity matrix; the vector of externally 

4> 



applied force; and F? the vector of externally applied charge, respectively 
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F « = lX N « p *i j i ded?? 



(29) 



(30) 



(31) 



(32) 



(33) 



(34) 



where p is the mass density, P s is the surface force vector, q s is the surface charge density vector. J 
is the well-known Jocobian matrix of the mapping (13) which is defined as: 

dx dy 



J = 



d(x,y) 



dx dy 
drj drj 



(35) 



The GLL integration rule is then used to calculate the characteristic matrices and the nodal force 
vector in Equation (28) at the elemental level [36]. In this study, the interface between the two PZT 
layers is grounded. Two sets of electric boundary conditions are considered, i.e., (1) sensor function 
with the top and bottom surfaces grounded and a uniform pressure load of S 0 = l,000N/m 2 applied 
to the upper surface, and Equation (2) actuator function with an electric potential of V Q = 50V 
applied to the top and bottom surface of the bimorph. By applying the electric boundary conditions, the 
DOFs for the electric potential are condensed out such that Equation (28) is finally of the form: 



uvr*- 1 uu 1 p ^ 



u 1 a 



(36) 



where denotes the mechanical stiffness matrix induced by the electromechanical coupling of PZT-4 
layer, and denotes the mechanical forces induced by the applied voltages of piezoelectric actuators [2]. 
The electric potential is then recovered by the inverse process of the aforementioned condensation. 
Assembling all elementary equations, one can have a global dynamic system equation: 



M ?m q+ K + K q 



F +F 

u 1 a 



(37) 



UU 



K lu 



where ]V[ „ ; , K „, , K , F , and F^ are the assembled counterparts of matrices M 

UU UU p U CL 

and F^ ; q is the global nodal displacement vector. Since the DOFs for the sublayer electric potentials 
have been condensed out, this approach will not result in a large number of potential DOFs. For the 
purpose of static analysis, the governing equations in Equation (37) reduces to: 
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K uu + K i 



F +F 

u ' a 



(38) 



3. Numerical Results 



In this section, the derived SE model is converted into a numerical code and case studies are carried 
out to validate the effectiveness and the capability of the present model for predicting both the global 
responses and the local responses, i.e., the deflections of the bimorph and the distribution of the 
electric potential across the bimorph thickness. A simply supported rectangular piezoelectric bimorph 
shown in Figure 1, which has been investigated by Fernandes [1], is considered here. The material 
constants of PZT-4 are given as: 
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(39) 
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(40) 



g 



13.06 0 0 
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nF/m 



(41) 



The length a and width b of the bimorph are 25 mm and 12.5 mm respectively. Two values of 
slenderness ratio, S = a / h = 5 and S = 50 , which represent the thick and thin bimorph plate, 
respectively, are considered. Unless otherwise stated, the order of Legendre polynomial is chosen as 5, 
and the mesh in Figure 3 is used in this work. Two load cases corresponding respectively to sensor 
function and actuator function are considered. To overcome the ill condition problem resulted from the 
huge difference of the element values of K e uu and in magnitudes, Equation (28) is rewritten using 
dimensionless variables. Consequently, the numerical results for the deflection and the electric 
potential are given in dimensionless units as: 



W,^ =^(w^/E Q ) for sensor function 



hS, 



E 



W,<& = —(w,4> I E 0 ) for actuator function 



(42) 



(43) 



where the amplification factor E Q is taken as E Q — 10 10 V/m . For the purpose of comparison, a 
coupled 3-D analysis is carried out using 20-noded hexahedral 3-D piezoelectric elements (C3D20RE) 
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with a mesh size of 40 x 20 x 10 in ABAQUS and the results from the coupled 3-D FE analysis are 
taken as accurate. 

3.1. Sensor Function 

For this case a uniform pressure load of S 0 = l,000N/m 2 is applied to the upper surface and the 
bimorph is used as a sensor with the top and bottom surfaces grounded. The variations of both the 
deflection W and the electric potential $ across the bimorph thickness at the centre of the bimorph 
plate (x = 0.5a, y = 0.56 ) for the slenderness ratio S = 5 and S = 50 are shown in Figures 5 and 6, 
respectively. It can be observed from Figure 5a and Figure 6a that the deflection W estimated by the 
present method adopting different number of sublayers is constant through the thickness and it is a 
good approximation of the nonlinear distribution described by the coupled 3-D analysis. The present 
model based on FSDT with assumption of uniform deflection through the thickness cannot predict the 
nonlinear variation of W through the thickness. The electric potentials induced by the deformation of 
the bimorph through the direct piezoelectric effects are shown in Figures 5b and 6b. It is observed that 
the distribution of the electric potential $ across the thickness provided by the present approach with 
more than 2 sublayers is in good agreement with the nonlinear distribution predicted by the coupled 
3-D analysis. Furthermore, it is expected that with more sublayers adopted the quadratic distribution of 
the electric potential $ across the bimorph thickness can be accurately approached without 
introducing any higher-order electric potential assumptions. However, the conventional linear electric 
potential assumption [38] will result in an inaccurate prediction of the local electric potential response 
for the case of sensors. The curves in Figures 5 and 6 are symmetrical with respect to the interface 
between the two PZT layers. It should be highlighted that although the present method cannot predict 
accurately the distribution of W across the bimorph thickness, it may be able to provide good 
approximate results for $ with appropriate number of sublayers for both thick and thin bimorph plates. 

Figure 5. Bimorph sensor of S = 5 under pressure load, (a) Dimensionless deflection; 
(b) Dimensionless electric potential. 3-D FE analysis (full line), present model with n = 2 
(triangles) and present model with n = 10 (small circles). 
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Figure 6. Bimorph sensor of S = 50 under pressure load, (a) Dimensionless deflection; 
(b) Dimensionless electric potential. 3-D FE analysis (full line), present model with n = 2 
(triangles) and present model with n = 10 (small circles). 





3.2. Actuator Function 

To achieve practically meaningful actuation capabilities and guarantee that the piezoelectric 
material behaves linearly, an electric potential of Vb = 50V is applied to the top and bottom surfaces of 
the bimorph with intermediate electrode grounded. The through-the-thickness variations of W and $ 
at the centre of the plate for S = 5 and S = 50 are shown in Figures 7 and 8, respectively. 

Figure 7. Bimorph actuator of S = 5 under potential load, (a) Dimensionless deflection; 
(b) Dimensionless electric potential. 3-D FE analysis (full line), present model with n = 2 
(triangles) and present model with n = 10 (small circles). 
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Figure 8. Bimorph actuator of S = 50 under potential load, (a) Dimensionless deflection; 
(b) Dimensionless electric potential. 3-D FE analysis (full line), present model with n = 2 
(triangles) and present model with n = 10 (small circles). 




(a) (b) 

Once again, the present model based on FSDT cannot predict accurately the through-the-thickness 
distribution of W. Similar to the previous observation, the constant deflection W through the thickness 
calculated by the present method adopting different number of sublayers is a good approximation of 
the nonlinear distribution provided by the coupled 3-D analysis. It is noticed that as the sublayer 
number increases a smaller deflection is obtained which is also pointed out by Wang [2]. The electric 
potentials at the centre of the plate are plotted in Figure 7b and Figure 8b for the slenderness ratio 
g = 5 and S = 50, respectively. It can be observed that the almost linear distribution of $ across the 
thickness predicted by the present method for both thick and thin bimorph plates is in excellent 
agreement with the coupled 3-D analysis, indicating that the nonlinear induced electric potential is 
insignificant compared to the externally applied potential. Consequently, the conventional linear 
electric potential assumption [38] may be accurate enough to calculate the local electric potential 
response for the case of actuators. 

4. Conclusions 

The present work aims to develop an efficient SE model with electric potential DOFs for the static 
electromechanical response of a piezoelectric bimorph. The approach is the combination of an ESL 
model based on FSDT for the mechanical displacement with a sublayer model based on the piecewise 
linear approximation for the electric potential. 2-D GLL shape functions are used to discretize the 
displacements and then the governing equation of motion is derived following the standard SEM 
procedure. By applying the electric boundary conditions, the DOFs for the electric potential are 
condensed out such that the present model will not result in a large number of potential DOFs. 

Numerical simulations based on the present model are carried out for two different load cases, i.e., a 
uniform pressure load applied to the top surface and a uniform potential applied to the top and bottom 
surfaces. To validate the effectiveness and the capability of the present model for investigation of both 
global and local response of the piezoelectric bimorph, the numerical results thus obtained are 
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compared to those from 3-D analysis using ABAQUS. The results indicate that the deflection W 
estimated by the present method is a good approximation of the nonlinear distribution predicted by the 
coupled 3-D analysis. It is further shown that the present model provides very accurate prediction for 
the electric potential distributions across the bimorph thickness even for rather thick bimorph plate 
without introducing any higher-order electric potential assumptions. It is also revealed that the 
conventional linear electric potential model is accurate enough to predict the local electric potential 
response for the case of actuators. This observation consists with the previous findings proposed by 
Yang [39]. One of the limitations is that the deflection W across the thickness is constant. 
Nevertheless, it is accurate enough to investigate the global response of the piezoelectric bimorph. The 
present work is important for researchers to better understand the nonlinear induced electric potential 
for bimorph sensor and actuator. An important extension of the present research is to study the 
vibration characteristics of the piezoelectric bimorph based on SE method. The influence of the 
induced stiffness matrix on the natural frequencies of the bimorph plate under various electric 
boundary conditions is to be investigated. The convergence study of the present model with respect to 
the order of the Legendre polynomial is also a practical and interesting problem to be conducted. 
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